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Abstract 

The distribution of resistance fluctuations of conducting thin films with granular 
structure near electrical breakdown is studied by numerical simulations. The film 
is modeled as a resistor network in a steady state determined by the competition 
between two biased processes, breaking and recovery. Systems of different sizes and 
with different levels of internal disorder are considered. Sharp deviations from a 
Gaussian distribution are found near breakdown and the effect increases with the 
degree of internal disorder. However, we show that in general this non-Gaussianity 
is related to the finite size of the system and vanishes in the large size limit. Never- 
theless, near the critical point of the conductor-insulator transition, deviations from 
Gaussianity persist when the size is increased and the distribution of resistance fluc- 
tuations is well fitted by the universal Bramwell-Holdsworth-Pinton distribution. 
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1 Introduction and Model 

Non-Gaussian distributions of several quantities characterizing the behavior of 
complex systems in non-equilibrium steady-states, have been evidenced since 
long time in many experiments [1]. The general properties of these anomalous 
distributions and their link with other properties of non-equilibrium complex 
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systems is far from being fully understood. Rather, the study of these distribu- 
tions is nowadays attracting an increasing interest in the literature [2,3,4,5,6,7]. 
Here, we study the distribution of resistance fluctuations of a conducting thin 
film with granular structure near electrical breakdown. Electrical breakdown 
phenomena typically occur in conductors stressed by high current densities 
and they consist in an irreversible and dramatic increase of resistivity of the 
material and thus they are associated with a conductor-insulator transition 
[8,9,10,11,12]. In our study we make use of the Stationary Network Under 
Biased Percolation (SNUBP) model [13]. This model provides a good descrip- 
tion of many features associated with the electrical instability of composites 
materials [11] and with the electromigration damage of metal lines [14], two 
important classes of breakdown phenomena. The film is modeled as a resistor 
network which reaches a steady state determined by the competition between 
two biased stochastic processes, breaking and recovery. Systems of different 
sizes and with different levels of internal disorder are considered, where the dis- 
order is mainly related to the fraction of broken resistors within the network. 
The resistance and its fluctuations are calculated by Monte Carlo simulations 
which are performed under different stress conditions. Resistance fluctuations 
are found to deviate from Gaussianity near electrical breakdown. We analyze 
and discuss this non-Gaussianity in the spirit of the Bramwell-Holdsworth- 
Pinton (BHP) distribution, recently introduced in the context of the study 
of highly correlated systems near criticality [3,4]. In particular, we show that 
the deviations from Gaussianity observed near the breakdown are related to 
finite size effects, and consequently vanish in the thermodynamic limit. Nev- 
ertheless, near the percolation critical point the non-Gaussianity persists in 
the large size limit and it is described by the universal BHP distribution. 
According to the SNUBP model [13], a conducting film with granular struc- 
ture is described as a two-dimensional resistor network. Precisely, we con- 
sider a square-lattice of N x N resistors. This lattice lies on an insulating 
substrate at a given temperature T , which acts as a thermal bath. Each 
resistor can be in two different states: (i) regular, corresponding to a re- 
sistance r n = ro[l + cx{T n — To)] and (ii) broken, corresponding to the ef- 
fectively "infinite" resistance, top = 10 9 r n (in the following we call this 
state defect). Here, a is the temperature coefficient of the resistance and 
T n the local temperature. This latter is determined by taking into account 
Joule heating effects and thermal exchanges between neighbor resistors [12]: 
T n = T + A[rJl + (3/4N neig ) E(n^ 2 - r n i 2 n )\, where i n is the current flow- 
ing in the nth resistor, N ne i g the number of nearest neighbors over which the 
summation is performed. The parameter A represents the thermal resistance 
of each resistor and determines the importance of Joule heating effects. By 
taking the above expression for T n we are neglecting time dependent effects in 
heat diffusion studied by Sornette et al. [9]. A constant stress current / is then 
applied through perfectly conducting bars at the left and right sides of the net- 
work. The two biased processes consist of stochastic transitions between the 
two possible states of each resistor and they occur, through thermal activa- 
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tion, with probabilities: Won = exp[— Eu/ksT^ and Wr u = exp[—EFt/kBT n ], 
characterized by the two activation energies, Ed and Er [13,15] (ks is the 
Boltzmann constant). The network time evolution is obtained by a Monte 
Carlo simulation which updates the network resistance after breaking and re- 
covery processes, according to the iterative procedure described in details in 
Ref. [13]. The sequence of successive configurations provides a resistance R(t) 
signal after an appropriate calibration of the time scale. Then, depending on 
the stress conditions (I and T ) and on the network parameters (sizes, acti- 
vation energies and other parameters related to the material like r and a), 
the network reaches a steady state or undergoes an irreversible electrical fail- 
ure. This last possibility is associated with the achievement of the percolation 
threshold, p c , for the fraction of broken resistors. Therefore, for a given net- 
work at a given temperature, a threshold current value, Ib, exists above which 
electrical breakdown occurs [13]. Below this threshold, the steady state of the 
network is characterized by fluctuations of the fraction of broken resistors, Sp, 
and of the resistance, SR, around their respective average values < p > and 
< R >. In particular, we underline that the ratio (Ed — Er) /ksT^ determines 
the average fraction of defects for a given value of I < Ib and thus the level 
of disorder inside the network. In the following section we analyze the results 
of simulations carried out by considering networks of different sizes, with dif- 
ferent levels of disorder and stressed by different currents. In all cases we take 
T = 300 (K), E D = 0.170 (eV), r = 1 (SI), a = 10~ 3 (K" 1 ), A = 5 x 10 5 
(K/W) (these values are chosen as physically plausible). N ranges between 
50 -r- 125, while E R between 0.0259 -=- 0.164 (eV). 



2 Results 

We report in Fig. 1 the resistance of a 75 x 75 network as a function of time and 
at increasing currents. The evolutions R(t) in Fig. 1 are obtained by taking 
Er = 0.103 (eV) for the activation energy of the recovery process. This value 
of Er provides a network with an intermediate level of disorder. This figure 
displays two important features of the electrical response of a conducting film: 
i) the linear regime at low currents is followed by a nonlinear regime where the 
average resistance increases significantly with the current; ii) the amplitude of 
the resistance fluctuations increases strongly close to the breakdown. Precisely, 
the lowest curve in Fig. 1 corresponds to the linear regime, i.e. it is obtained 
for I = 0.001(A) < Iq where Jo is the current value associated with the onset 
of the nonlinearity of the I-V characteristic [13]. The second curve corresponds 
to the nonlinear regime, I = 0.45(A) > I , the third to the threshold for the 
breakdown, I — I B — 0.70(A), and finally the highest curve to a network 
undergoing an electrical breakdown at a current I = 0.75(A) > Ib- This last 
curve has been shifted upwards by 0.1(1)) for graphical reasons. A detailed 
analysis of the behavior of the average resistance and of the relative variance 
of resistance fluctuations as a function of the current can be found in Ref. [13]. 
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Fig. 1. Resistance evolutions at increasing currents. Starting from the bottom: 
/ = 0.001 (A) (linear regime), / = 0.45 A) (nonlinear regime), / = 0.70 (A) (thresh- 
old current), / = 0.75 (A) (breakdown). The highest curve has been shifted upwards 
by 0.1(0) for graphical reasons. 
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Fig. 2. Normalized PDF of resistance fluctuations in the linear regime at increasing 
network sizes: 50 x 50 (stars), 75 x 75 (small circles), 100 x 100 (squares), 125 x 125 
(down triangles). The recovery energy is En = 0.103 (eV). The thick solid curve and 
the dashed curve correspond to the BHP and Gaussian distributions, respectively. 



We have investigated the Gaussianity of the steady state signals shown in Fig. 
1. The results of the analysis are reported in Fig. 2 for / = 0.001 (A) (linear 
regime) and in Fig. 3 for I = I B = 0.70 (A) (threshold current for breakdown). 
Precisely, in these figures we show on a lin-log plot the product <r$ as a function 
of (< R > —R)/a, where $ is the probability density function (PDF) of the 
distribution of SR and a is the root mean square deviation from the average. 
This normalized representation, by making the distribution independent of its 
first and second moments, is particularly convenient to explore the functional 
form of the distribution. These results, together with all the others presented 
in this paper have been obtained by considering time series containing about 
1.2 x 10 6 resistance values. Moreover, for comparison, we have also reported the 
Gaussian distribution (which in this normalized representation has zero mean 
and unit variance) and the Bramwell-Holdsworth-Pinton (BHP) distribution 
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Fig. 3. Normalized PDF of resistance fluctuations at the threshold current, Ib, at 
increasing network sizes: 50 x 50 (stars), 75 x 75 (small circles), 100 x 100 (squares), 
125 x 125 (down triangles). The recovery energy Er is the same of Fig. 2. The thick 
solid and the dashed curves also have the same meaning. 

[3,4]. These authors have proposed a universal non-Gaussian PDF for the 
fluctuations of global quantities of systems at criticality[3]. Denoting to a 
fluctuating quantity (for example the magnetization of a ferromagnet like in 
Ref. [4]), < to > and a m its mean value and root mean square deviation 
respectively, P(m) its PDF, y = (to— < to >)/cr m the normalized variable, 
Il(y) = a m P(y) the normalized PDF and x = b(y — s), the BHP distribution 
takes the following expression [4]: 



where a = vr/2, b = 0.936±0.002, s = 0.374±0.001 and K = 2.15±0.01. This 
expression can be considered as a generalization of the Gumbel distribution, 
which is often associated with the occurrence of rare events. In Figs. 2 and 
3 together with the PDF obtained for networks with linear sizes N = 75 
(small circles) are also displayed the PDF corresponding to iV = 50 (stars), 
N = 100 (squares) and N = 125 (down triangles). Figure 2 shows that in 
the linear regime the distribution of SR is Gaussian for all system sizes. Near 
breakdown, at I = Ib, (Fig. 3) the PDF shows non-Gaussian tails. However, 
this non-Gaussianity progressively vanishes when the system size is increased. 
Therefore, when networks with intermediate level of disorder are considered, 
as in the case of Figs. 2 and 3, deviations from Gaussianity are related to the 
finite size of the system. The PDF in this case is well fitted by the Eq. (1) 
once the parameters a, b, s and K are considered as fitting parameters, as 
described in Ref. [16]. Such phenomenological fits can provide helpful tools in 
the study of failure precursors in the case of finite size systems. 
We have also considered the role of disorder in the breakdown process and its 
effect on the distribution of resistance fluctuations. At given values of E D and 
To, the parameter controlling the disorder inside the network, i.e. the average 
fraction of broken resistors < p >, is the activation energy Er. Figure 4 shows 
the average fraction of broken resistors versus the external current for different 
values of E R . The data are obtained for a 75 x 75 network. The different sets of 
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Fig. 4. Average fraction of defects versus current in steady-state networks for dif- 
ferent values of the recovery energy Er. The numbers at the top of each curves 
indicate Er in eV. In the region above the solid curve the system breaks down. 
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Fig. 5. Normalized PDF of resistance fluctuations of a 75 x 75 network at increasing 
Er values. In all cases the current corresponds to the threshold value. Precisely, left 
triangles: Er = 0.026, Ir = 2.6; full triangles: Er = 0.043, Ib = 2.1; small circles: 
E R = 0.103, I B = 0.70; diamonds: E R = 0.140, I B = 0.25; full squares: E R = 0.155, 
Ib = 0.085, crosses: Er = 0.164, Ir = 0.0065 (energies in eV, currents in A). 



symbols refer to steady states of the network. The region above the solid curve 
corresponds to electrical breakdown (in this region the stable states of the 
network are non-conducting states). Thus, Fig. 4 represents the phase diagram 
of the system in the < p > —I plane. It can be shown [17] that all the curves 
in Fig. 4 collapse onto a single one after normalization of / to 1$ and of < p > 
to < p > (the defect fraction in the limit of vanishing current). Moreover, the 
relative variation of the defect fraction, [< p > — < p > ]/ < p > , scales as 
the ratio I/I with a quadratic exponent [13]. Furthermore, we note that the 
electrical breakdown corresponding to the values of E R considered in Fig. 4, is 
associated with a first order phase transition [17]. These results agree with the 
behavior observed in electrical breakdown experiments, performed in the Joule 
regime of composites [11]. Nevertheless, it can be shown [17] that when Er 
reaches its maximum value, E RtM Ax, (determined by the stability condition 
obtained for a system of given size in the vanishing current limit [15,17]), 
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Fig. 6. Normalized PDF of resistance fluctuations at the threshold current, Ir, at 
increasing network sizes: 75 x 75 (crosses), 100 x 100 (down triangles), 125 x 125 
(full circles). The recovery energy is Er = 0.164 (eV). 

the conductor-insulator transition becomes of the second order. This change 
in the nature of the transition, when going from small to high disordered 
systems, also agrees with the predictions of Andersen et al. [9]. Therefore, 
we have reported in Fig. 5 the normalized PDF calculated at the breakdown 
threshold for different values of Er. Again, a 75 x 75 network is considered. 
We can see that the non-Gaussianity of the distribution at I = Ib increases 
significantly when increasing the Er value. However, as already shown in Fig. 
3, this non-Gaussianity vanishes in the large size limit. Nevertheless, when 
the value of E R is nearly equal to E r>M ax, the system approaches the critical 
point and the PDF tends to achieve the BHP form, independently of system 
size. This is shown in Fig. 6 which displays the normalized PDF calculated at 
the breakdown threshold, I = Ib, for Er = 0.164 (eV). The three sets of data 
are obtained for N = 75, 1 = 0.0065 (A) (crosses), N = 100 and I = 0.009 (A) 
(down triangles) and iV = 125 and I = 0.011 (A) (full circles) and correspond 
to networks near their critical point (the data for N = 50 are not shown in 
Fig. 6 because 50 x 50 networks are always unstable for this value of Er). 
We can see that in this case the PDF becomes independent of the sizes of the 
system and it is well fitted by the BHP distribution. 

3 Conclusions 

We have studied the distribution of the resistance fluctuations of conducting 
thin films with granular structure near electrical breakdown. The study has 
been performed by describing the film as a resistor network and by using 
the SNUBP model [13]. We have considered systems of different sizes and 
with different levels of internal disorder. A non-Gaussianity of the fluctuation 
distribution is found near electrical breakdown. This non-Gaussianity increases 
with the degree of disorder of the network. However, we have shown that this 
non-Gaussianity is related to the finite size of the system and that it vanishes 
in the large size limit. Nevertheless, near the critical point the non-Gaussianity 
persists in the large size limit and is well fitted by the universal Bramwell- 
Holdsworth-Pinton distribution [3,4,5]. 
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